from numpy import pi, zeros, sin, cos, arctan, linspace
from bem import TwoDimensionalLaplaceBEMSolver

NO = 5

N = 12 * NO

xb = zeros(N + 1)
yb = zeros(N + 1)

bct = zeros(N)
bcv = zeros(N)

for i in range(8 * NO):

    dL = pi / float(16.0 * NO)

    theta = float(i) * dL

    xb[i + NO] = 2.0 * cos(theta)
    yb[i + NO] = 2.0 * sin(theta)

    if i < NO:
        dL = 1.0 / float(NO)
        xb[i] = 1.0 + float(i) * dL
        yb[i] = 0.0
        xb[i + 9 * NO] = 0.0
        yb[i + 9 * NO] = 2.0 - float(i) * dL

    if i < 2 * NO:
        dL = pi / float(4 * NO)
        theta = float(i) * dL

        xb[i + 10 * NO] = sin(theta)
        yb[i + 10 * NO] = cos(theta)

xb[N] = xb[0]
yb[N] = yb[0]

bem = TwoDimensionalLaplaceBEMSolver(xb, yb)

for i in range(N):
    if i < NO or i >= 9 * NO and i < 10 * NO:
        bct[i] = 1
        bcv[i] = 0.0
    elif i >= NO and i < 9 * NO:
        bct[i] = 0
        bcv[i] = 3.0 * cos(4.0 * arctan(bem.ym[i] / bem.xm[i]))
    else:
        bct[i] = 0
        bcv[i] = cos(4.0 * arctan(bem.ym[i] / bem.xm[i]))

bem.bctypes = bct
bem.bcvalues = bcv

bem.solve()

points = [(1.082532, 0.625000),
          (0.875, 1.515544),
          (1.06066, 1.06066),
          (1.09998, 0.00192),
          (1.01, 0.000176)]

for point in points:
    print bem.evaluate_phi_at_interior_point(point[0], point[1])


import matplotlib.pyplot as plt

npoints = 100
r = linspace(1.1, 1.9, npoints)
theta = linspace(0.0, pi / 2.0, npoints)

rr = linspace(1.0, 2.0, npoints)

x1 =  cos(theta)
y1 =  sin(theta)

x2 = 2.0 * cos(theta)
y2 = 2.0 * sin(theta)

xx = zeros((npoints, npoints), float)
yy = zeros((npoints, npoints), float)
zz = zeros((npoints, npoints), float)

for m in range(r.size):
    for n in range(theta.size):
        xx[m, n] = r[m] * cos(theta[n])
        yy[m, n] = r[m] * sin(theta[n])
        zz[m, n] = bem.evaluate_phi_at_interior_point(xx[m,n], yy[m,n])

plt.plot(x1, y1, 'black')
plt.plot(x2, y2, 'black')

plt.contourf(xx, yy, zz)
plt.grid()
plt.show()
